{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Install Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: numpy in /home/rzhang/anaconda3/lib/python3.6/site-packages (1.16.1)\n",
      "Requirement already satisfied: scipy==1.2.0 in /home/rzhang/anaconda3/lib/python3.6/site-packages (1.2.0)\n",
      "Requirement already satisfied: matplotlib in /home/rzhang/anaconda3/lib/python3.6/site-packages (3.1.0)\n",
      "Requirement already satisfied: autograd==1.1.13 in /home/rzhang/anaconda3/lib/python3.6/site-packages (1.1.13)\n",
      "Requirement already satisfied: tick==0.5.0.0 in /home/rzhang/anaconda3/lib/python3.6/site-packages (0.5.0.0)\n",
      "Requirement already satisfied: numpydoc==0.7.0 in /home/rzhang/anaconda3/lib/python3.6/site-packages (0.7.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from matplotlib) (2.4.0)\n",
      "Requirement already satisfied: cycler>=0.10 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from matplotlib) (0.10.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from matplotlib) (1.1.0)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from matplotlib) (2.8.0)\n",
      "Requirement already satisfied: future>=0.15.2 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from autograd==1.1.13) (0.17.1)\n",
      "Requirement already satisfied: sphinx in /home/rzhang/anaconda3/lib/python3.6/site-packages (from tick==0.5.0.0) (1.7.4)\n",
      "Requirement already satisfied: pandas in /home/rzhang/anaconda3/lib/python3.6/site-packages (from tick==0.5.0.0) (0.23.0)\n",
      "Requirement already satisfied: scikit-learn in /home/rzhang/anaconda3/lib/python3.6/site-packages (from tick==0.5.0.0) (0.19.1)\n",
      "Requirement already satisfied: Jinja2>=2.3 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from numpydoc==0.7.0) (2.10)\n",
      "Requirement already satisfied: six in /home/rzhang/anaconda3/lib/python3.6/site-packages (from cycler>=0.10->matplotlib) (1.12.0)\n",
      "Requirement already satisfied: setuptools in /home/rzhang/anaconda3/lib/python3.6/site-packages (from kiwisolver>=1.0.1->matplotlib) (41.0.1)\n",
      "Requirement already satisfied: Pygments>=2.0 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (2.2.0)\n",
      "Requirement already satisfied: docutils>=0.11 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (0.14)\n",
      "Requirement already satisfied: snowballstemmer>=1.1 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (1.2.1)\n",
      "Requirement already satisfied: babel!=2.0,>=1.3 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (2.5.3)\n",
      "Requirement already satisfied: alabaster<0.8,>=0.7 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (0.7.10)\n",
      "Requirement already satisfied: imagesize in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (1.0.0)\n",
      "Requirement already satisfied: requests>=2.0.0 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (2.18.4)\n",
      "Requirement already satisfied: packaging in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (17.1)\n",
      "Requirement already satisfied: sphinxcontrib-websupport in /home/rzhang/anaconda3/lib/python3.6/site-packages (from sphinx->tick==0.5.0.0) (1.0.1)\n",
      "Requirement already satisfied: pytz>=2011k in /home/rzhang/anaconda3/lib/python3.6/site-packages (from pandas->tick==0.5.0.0) (2018.4)\n",
      "Requirement already satisfied: MarkupSafe>=0.23 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from Jinja2>=2.3->numpydoc==0.7.0) (1.0)\n",
      "Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from requests>=2.0.0->sphinx->tick==0.5.0.0) (3.0.4)\n",
      "Requirement already satisfied: idna<2.7,>=2.5 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from requests>=2.0.0->sphinx->tick==0.5.0.0) (2.6)\n",
      "Requirement already satisfied: urllib3<1.23,>=1.21.1 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from requests>=2.0.0->sphinx->tick==0.5.0.0) (1.22)\n",
      "Requirement already satisfied: certifi>=2017.4.17 in /home/rzhang/anaconda3/lib/python3.6/site-packages (from requests>=2.0.0->sphinx->tick==0.5.0.0) (2018.4.16)\n",
      "\u001b[33mYou are using pip version 18.0, however version 19.3.1 is available.\n",
      "You should consider upgrading via the 'pip install --upgrade pip' command.\u001b[0m\n",
      "Cloning into 'VBHP'...\n",
      "remote: Enumerating objects: 31, done.\u001b[K\n",
      "remote: Counting objects: 100% (31/31), done.\u001b[K\n",
      "remote: Compressing objects: 100% (30/30), done.\u001b[K\n",
      "remote: Total 31 (delta 9), reused 0 (delta 0), pack-reused 0\u001b[K\n",
      "Unpacking objects: 100% (31/31), done.\n"
     ]
    }
   ],
   "source": [
    "!pip3 install numpy scipy==1.2.0 matplotlib autograd==1.1.13 tick==0.5.0.0 numpydoc==0.7.0\n",
    "import os \n",
    "if os.path.exists('VBHP'):\n",
    "    # update VBHP repo\n",
    "    ! rm -rf VBHP\n",
    "try:\n",
    "    !git clone https://github.com/RuiZhang2016/Variational-Inferencce-for-SGP-Modulated-Hawkes-Process.git VBHP\n",
    "except Exception as e:\n",
    "    print(e)\n",
    "import sys\n",
    "sys.path.append('./VBHP')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Import Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rzhang/anaconda3/lib/python3.6/site-packages/autograd/core.py:120: UserWarning: \n",
      "------------------------------\n",
      "    defgrad is deprecated!\n",
      "------------------------------\n",
      "Use defvjp instead (\"define vector-Jacobian product\").\n",
      "The interface is a little different - look at\n",
      "autograd/numpy/numpy_grads.py for examples.\n",
      "\n",
      "  warnings.warn(defgrad_deprecated)\n",
      "/home/rzhang/anaconda3/lib/python3.6/site-packages/autograd/core.py:120: UserWarning: \n",
      "------------------------------\n",
      "    defgrad is deprecated!\n",
      "------------------------------\n",
      "Use defvjp instead (\"define vector-Jacobian product\").\n",
      "The interface is a little different - look at\n",
      "autograd/numpy/numpy_grads.py for examples.\n",
      "\n",
      "  warnings.warn(defgrad_deprecated)\n",
      "/home/rzhang/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import autograd.scipy.special as special\n",
    "from util import *\n",
    "import equations\n",
    "import autograd.numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from autograd import grad\n",
    "from vbhp import fit_vbhp\n",
    "import wp\n",
    "from scipy.optimize import minimize\n",
    "from tick.hawkes import HawkesKernelTimeFunc, SimuHawkes\n",
    "from tick.base import TimeFunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define a Triggering Kernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a triggering kernel: $$\\phi(t)=\\begin{cases}0.9[\\sin(3t)+1], & t\\in [0,\\pi) \\\\ 0, & \\text{otherwise} \\end{cases}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = np.linspace(0, np.pi / 2, 512) # support\n",
    "f = lambda t: (np.sin(3 * t)+1)*0.9\n",
    "ys = f(ts)\n",
    "kernel = HawkesKernelTimeFunc(t_values=ts, y_values=ys) # for simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### plot the triggering kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEaCAYAAACfC2mcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZgU1bn48e87G9uwCMyw7yACgoADbhFU3EVFRU3czWKue35ZXGJurjfJFdSbqInRxKsGjUtUXAm4gIK7yLDvAsq+DTvDDvP+/jg1M9093T09M71UT7+f5+mnqFOnqt/umuGdOnXqHFFVjDHGmHSUleoAjDHGmNqyJGaMMSZtWRIzxhiTtiyJGWOMSVuWxIwxxqStnFQHkGlat26tXbt2TXUYxhiTNmbOnLlFVQvCbbMklmRdu3aluLg41WEYY0zaEJFVkbZZc6Ixxpi0ZUnMGGNM2rIkZowxJm1ZEjPGGJO2LIkZY4xJW5bEjDHGpC3rYm8SrmT3AaYt3czyzaV8u2UPew8epkXjPFo1yWNAxxac3ruAVvkNUh2mMSYN+TaJiUhv4FxgCFAEHA0IcLmqjq/hsU4DpsZYvYuqrg7YdxxwfZT6S1X1mJrEkwlUlWlLS3hx+iqmLi3hSFmkKX9WkSUwpGtLbjm9J8N6tUZEkhqrMSZ9+TaJATcDd8bpWBuB56JsHwr0AVYAayLU+RxYHqZ8Q91Cq39WlJRy/zsL+XTZlpjqlylM/24b07/7mhO6teQ3F/Slf8fmCY7SGFMf+DmJLQAeBoqBmcAzwPDaHEhVlwA3RNouIou8fz6rkWcJfVpVx9Xm/TOFqvJ/n37Lw+8v5dCR2k22Ov27bVzyxOfcc94x/Oh73eyqzBgTlW+TmKo+HbieqP/MROQk3FXYEWBcQt4kA+w/dIR7Xp/HW3PWh91+dJt8zurbhp6F+RzVOI+d+w6xfHMpHy7ezKINu4LqHi5T/jBxMV99u5VHrhxI04a5yfgIxpg05NsklkQ/9JbvqWr4/4FNVDv3HuL6f3zNnDU7qmw7qXsr7jq3NwM7tQj7h8gvzu7N/LU7efiDpXzyTUnQtimLN3PNM1/z/I1Dad7YEpkxpqqMTmIi0hi40lt9pprqp4vIACAf2AR8BkxW1bIEhuh7u/cf4rp/fM3ckATWrGEOf7ikPxcOaFftVXT/js15/odDeX/hRn712lx27T9csW3umh384P++4p8/Gmo9GI0xVWR0EgMuB5oCm4F/V1P3ujBli0Tk+6o6P+6RpYG9Bw/zo3HFVRJYj4ImPH39ELq1blKj453Try192zXjlhdnMX/dzoryRRt2ce0zX/Paf5xEkwaZ/iNrjAmU6Q87lzclPq+qhyLUmQPcAfTFXYW1B0YCc72yKSLSIdqbiMhNIlIsIsUlJSXRqqaNsjLlzn/N4euV24LKh3ZtyZu3nlLjBFauU8vGvHzTiQzt1jKofNGGXdz5r9lRuuobYzJRxiYxEekJDPNWn41UT1UfVdW/qOpiVd2jqhtUdSKuW/5XQCFwb7T3UtWnVLVIVYsKCsLO65Z2/vzRMiYv2hRUNrBTC569cQjN6tgRI79BDs/dOJTv9WwdVD5l8WbGTFpcp2MbY+qXjE1iVF6FfamqNf6fUVUPAmO81fPjFlUamLxoE49OWRZU1qddM5774VDy49Tc1ygvm79fezx92zULKn/6s+94e866uLyHMSb9ZWQSE5FsKu9xVdehI5ol3jJqc2J9snb7Xn7+ypygslZN8njm+iKaN4pvD8ImDXJ45oYi2jQL7tBx35sLWL11b1zfyxiTnjIyiQHn4BJPKfBKHY7TyluW1jmiNFBWpvzytbnsPlDZezA7S/jr1YNp36JRQt6zXfNGPH3dEPJyKn9USw8c5s5XZnPoSEZ3DDXGkLlJ7Efe8lVVrUsCusJbzqhjPGnh2c+/46tvgzty/Pr8PpzYvVWEPeKjf8fm3Hd+n6Cy2at38JcPl0XYwxiTKepVEhORMSKyRETGRKnTGrjQW43alCgiA0VkpNf8GFieIyK/wPVaBHikLnGng+Wbd/PQ+0uDyk7t1ZofntI1Ke9/3UldGHFMYVDZE9NWsHTj7qS8vzHGn3ybxERksIh8Vf4CBnubHggpD9QO6O0tI7kWyAWWqOoX1YTRFZgAbBaRySLyooi8B6wC/terc5eqvl+Dj5Z2VJVfv7GAg4crm++aNczhodEDkja2oYjw0OgBFDStvD92uEy55415lFm3e2Mylm+TGNAMOCHg1dQr7xVSXlM3esuI3eoDzAUeA5bingm7DDcI8V7gH8BQVX24FjGklTdmravyPNjvRx1Lu+aJuQ8WSav8Bvz+4n5BZbNX7+DF6auSGocxxj8k8qDtJhGKioq0uLg41WHEbOfeQ4z40zS2lB6sKDu9dwHP3jAkJSPMqyo3/XNm0DNq+Q1y+OiXwyls2jDp8RhjEk9EZqpqUbhtfr4SMz7wx8lLgxJYg5ws/vuiY1M2RYqI8N8X9aNJXuVtytIDh/nTB9+kJB5jTGpZEjMRLd+8mxe+Cm6qu+W0nnRu1ThFETntWzTiF2f3Dip7pXgNi9bvirCHMaa+siRmIhr77lIC+0x0adWYnw7vnrqAAlx7Uhe6F1SOz6gK/zNpEdY8bkxmsSRmwvr6u21MWRw8NuK95x1Dw9zsCHskV252VpVnxz5fvpWPlmxOUUTGmFSwJGaqUFUeCBlod3DnFpzTr22KIgrvjGMKOaVn8IPWD7+/1LrcG5NBLImZKt5fuKnKLM2/Pr9PyjpzRCIi3Hd+XwLDWrJxN5MWbEhdUMaYpLIkZoKUlSmPhQzndHbfNhR1bRlhj9Tq274ZFw5oH1T26JRlNu+YMRnCkpgJMnnxJhZvqOzlJwK/Oqd3lD1S784ze5EVcDW2fHMp78y16VqMyQSWxEwFVeXPIVdh5/dvR682TSPs4Q89CvK5ZFDHoLLHpizjsI1yb0y9Z0nMVJiyeDMLQ561uuOMXimKpmbuHNGLnIDLsZVb9/Lugo0pjMgYkwy+TWIi0ltE7hSRF7yR6ctEREVkdC2PN87bP9JrSZR9s0TkVhEpFpFSEdkpIp+KyA9q/wn9RVV5/KPQq7C29G7r76uwcp1bNebSwcFzkz4xbYU9N2ZMPRefueQT42bgzgQc93NgeZjysF3avGlY3gAuAnYBHwANgBHASyJyoqomIs6k+urbbcxduzOo7I4R6XEVVu6nw3vw2sy1lOetxRt28fE3JZzWuzD6jsaYtOXnJLYAeBgoBmbi5v4aHofjPq2q42pQ/2e4BLYIOENVNwGISC/gU+AOEflIVd+OQ2wp8/dPVgStn9mnkGPaNktRNLXToyCf845ty6T5lc2IT0xbYUnMmHrMt82Jqvq0qt6lqq+q6orq94g/7yrsLm/15vIE5sW3DLjbW70v2bHF05KNu5i2tCSo7KfDe6Qomrq5eXjPoPWvv9vGzFXbItQ2xqQ73yYxnzgJKATWquonYba/BhwChohIhzDb08JTn3wbtD64cwuKuhyVomjqpn/H5pzaq3VQ2TOffZeiaIwxiZaJSex0EfmTiDwlIr8XkXNEJNL3MMhbzgi3UVX3Agu91YHxDjQZNu7czztz1geV3TSsh+9G56iJnw4Lvop8b8FG1m7fm6JojDGJlIlJ7Drg/wE/AX4DvAfMF5H+Yep285bRpg5eHVI3rbzw1SoOB4xu0b11E87q2yaFEdXdKT1bcUxAr8oyhee/tNmfjamPMimJzQHuAPoC+UB7YCQw1yubEqZJMN9b7oly3FJvGbEvuojc5HXPLy4pKYlULen2HzrCS1+vDiq78XvdyM5K36swcGMq/vCU4L8pXv56NXsOHE5RRMaYRMmYJKaqj6rqX1R1saruUdUNqjoRGAp8hbv3dW+C3vspVS1S1aKCgoJEvEWtTJi7nm17Kmdtbtowh0sHpe2tvSAXDWxPyyZ5Feu79x/m9VlrUxiRMSYRMiaJRaKqB4Ex3ur5IZvLr7KaEFn51drueMaVaKrKuC9WBpVdWdSJJg38/NRF7BrmZnPNCZ2DysZ9vtIefjamnsn4JOYpH60j9DJkpbfsEmXfTiF100Lxqu1BQ0yJwHUndU1dQAlwzYldyM2ubBr9dssevlixNYURGWPizZKYUz6zYmlI+SxvOSTcTiLSGDjWW52dgLgS5p8hHR1GHNOGzq0apyiaxChs1pBzj20XVBb6uY0x6c2SmHOFtwztSv8lUAJ0FJFhYfa7HMgFZqhq2sz9sX3PQd4LGRz3hpO7piaYBLv2xOCL6MmLN7Fx5/4URWOMibd6lcREZIw3WPCYkPKBIjLSG4EjsDxHRH6B67UI8EjgdlU9AjzkrT4pIoUB+/YCxnqr/xPPz5Fob8xex8GAaUq6tGrMyT1aRdkjfQ3pehRHt8mvWD9Sprwc0iPTGJO+fJvERGSwiHxV/gIGe5seCCkP1A7o7S0DdQUmAJtFZLKIvCgi7+Ge//pfr85dqvp+mFAe8fbtCywTkTdEZAIwD2gL/CWdxk1UVf4V8p/4lUM6kZXm3eojEZEqV2Mvf72aQzbXmDH1gm+TGNAMOCHgVf4cVq+Q8ljMBR4DluKS0WW4wYT3Av8Ahqrqw+F29K7GRgG340a/P8fbdyZwtareEW4/v5q1ejvLNlfe+svJEkYf3zHKHulv1KAONM6rvAjfvPsAHy7enMKIjDHx4tv+1Ko6DajR5YGq3gDcEKb8O9xo9LWNpQx43HultZe/XhO0PqJPIYVNG6YomuRo2jCXUYM68NL0yivQ8TPXcO6xbVMYlTEmHvx8JWbibNf+Q/x7XvA4id8f0jlC7frlyqJOQetTl5awebd18DAm3VkSyyBvz1nP/kOV94LaN2/IsKP9M4JIIg3o2JzebSpHBjtSprw5K206lBpjIrAklkFCO3RcXtQp7cdJjJWIcHlR8L2/V4vX2AgexqQ5S2IZYv7anVVG6LhiSKcoe9Q/lwzqQE5A0l5RsodZq3ekMCJjTF1ZEssQL88IvgobfnQBHVo0SlE0qdEqvwEj+hQGlb1WvCZCbWNMOrAklgH2HDhcZeLLTOnQEeqKkA4eE+auZ+9Bm6LFmHRlSSwDTJq/gdKAubRah7kiyRTDjy6goGmDivU9B48waf7GKHsYY/zMklgGeHN2cC+80cd3JDc7M099TnYWlw0O7uBhTYrGpK/M/J8sg2zYuY8vvw2efmT08fVj4svaCu2lOP27bazcEm3ybmOMX1kSq+femr2ewF7k/Ts0p2dh08g7ZIAeBfkc3+WooLLxM23WZ2PSkW+TmIj0FpE7ReQFb2T6MhFRERldi2PlisgIEfmjiBSLyC4ROSgi60RkvIicFmXfcd77RnotibRvqqkqb84O/s/5kkGZfRVW7oqQq7HxM9dypMyeGTMm3fh27ETgZuDOOB1rODDZ+/dG4BNgD5WDAV8mIr9X1d9GOcbnuAGAQ22IU4xxt2jDLr7ZVDnYb3aWcNHA9imMyD8uGNCe+99ZxL5DRwDYuGs/07/dysk9W6c4MmNMTfg5iS0AHgaKcSPGP4NLRrVRBrwOPKaqnwZuEJErgReB/xSRqao6NcIxnlbVcbV8/5QIHVZpWK/WtM5vEKF2ZslvkMO5x7YN6vTy9pz1lsSMSTO+bU5U1adV9S5VfVVVV9TxWB+p6ujQBOZtewUY561eU5f38ZPDR8p4e27ws2GjrCkxyMUhV6WTFmxgv3dlZoxJD75NYkk221vWm4m1vlixlZLdByrW8xvkcHZfm3ok0Pd6tqZVk7yK9d37DzNtqc0zZkw6sSTm9PKW0e5vnS4ifxKRp0Tk9yJyjoj49vsLfTbs3GPb0ihgYkjjnhkbOSB4EvC3Zq+PUNsY40d+vieWFCLSlsqJNF+PUvW6MGWLROT7qjo/7oHVwZ4Dh3lvQfAoFJdaU2JYFw/qwHNfrqpY/2jpZnbuO0TzRrkpjMoYEyvfXkkkg4jkAC8AzYEPVXVCmGpzgDtwPRnzgfbASGCuVzZFRHyVId5fuLGi1x1A22YNOaF7qxRG5F+DOrWgS6vGFesHD5fx/gIbhsqYdJHRSQz4GzACWEOETh2q+qiq/kVVF6vqHlXdoKoTgaHAV0AhcG+0NxGRm7zn04pLSkri/BGqCm1KvHhQ+4yZN6ymRISLjwvu4PHWHJss05h0kbFJTEQeA36Ee25shKrW6M9vVT0IjPFWz6+m7lOqWqSqRQUFiZ1JefOu/Xy+fEtQ2aWD6k1/lYS4aGDwhfSX325l4879KYrGGFMTGZnEROSPuCbCElwCW1bLQ5WP1uGb5sSJ8zcQOPBEn3bN6N02s4eZqk7PwnyO7dCsYl3VTdFijPG/jEtiIvIQ8HNgK3Cmqi6qw+HKbzSVRq2VRP+eF9zBcpSN0BGTUSFXY2/PtSZFY9JBRiUxERkL/ArYDpylqvPqeMgrvOWMOh4nLtZu38vMVduDyi4I6UJuwrvwuPZIwG3DBet2sXyzb/42McZEUK+SmIiM8QYLHhNm2x+Au4EduAQ2u8oBqu4zUERGikh2SHmOiPwC1yQJ8Egcwq+ziSFXYYM7t6DjUY0j1DaB2jRryMk9gntwvm0dPIzxPd8+JyYig4EnAor6essHROSX5YWqemJAnXZAb28ZeKyLgPu81eXA7SJhe+stUdWxAetdgTeBbSIyC9iMa0Lsj+tqXwbcparv1+jDJUhoU+LIAdaUWBMXH9eBz5dXzr32ztz1/Pyso4nws2KM8QHfJjGgGXBCmPJeYcqq0zLg30XeK5yPgcAkNhd4DNedvi9wKqDAWuAfwF9VdWYt4om7lVv2MH/dzop1EWtKrKlz+7flN28v4ODhMgBWbd3LwvW7OLZD8xRHZoyJxLdJTFWnATX6E1hVb6By9I3A8nFUDvJbk+N9B/yspvulwr/nBfemG9q1JW2aNUxRNOmpWcNchh9dwORFmyrKJsxbb0nMGB+rV/fEMtmEucFNiRceZ02JtRE6luLEeRtQtckyjfErS2L1wDebdrN00+6K9ews4bxjbcT62hjRpw0Ncip/LdZu38e8tTuj7GGMSSVLYvXAv0MezD25Ryta2eSXtZLfIIfTexcGlU2c79vJu43JeJbE0pyqMmGeNSXGU2iHGGtSNMa/LImluYXrd/Hdlj0V67nZwjk2+WWdnHFMIQ1zK3811u3Yx+w1O1IYkTEmEktiaW5CSK/E4UcX0LyxzYVVF00a5DDimDZBZaEPkhtj/MGSWBpT1Sr/uVpTYnyEa1IsK7MmRWP8xpJYGpu3didrt++rWG+Qk8WIPm2i7GFidXrvQhrnVY42tnHXfmat3h5lD2NMKlgSS2OTFgRfhZ3eu5D8Br59fj2tNMrLrvIHQeiwXsaY1LMklqZUlUkhXb/P628dOuLpgv7BTYqT5luTojF+Y0ksTS1cv4s12yqbEvOsKTHuTutdQJOAJsXNuw8wY+W2FEZkjAnl2yQmIr1F5E4RecGbXqVMRFRERtfxuFeJyKcislNESkWkWERuFZGo34WInCsiH4jINhHZKyILROQ+EUnJU8WhV2HDjy6wpsQ4a5ibzVl9Q3op2oPPxvhKrf7XE5FWwOnAIKAN0AI30eRmYBYwTVW3Rj5CTG4G7qzjMYKIyF+BW4D9wIfAIWAE8DgwQkRGq2pZmP3uAh4EjgDTcJ91OPAHYKSIjFDVvfGMNZpwTYnnW1NiQlwwoD1vzal8jOHdBRv5rwv7kZ1l07MY4wcxJzERyQEuxyWBk3AjzIf7TVZAReQL3Hxg41X1cC1iWwA8DBQDM4FncImjVkTkMlzsG4FhqrrMK28DTAUuAW7HTb0SuF8RbnqWvcAZqjrdK88HJgLDgP8B/l9tY6upxRt2s3JrZc7My7amxEQ5tVdr8hvkUHrA/QiX7D5A8cptnNC9VTV7GmOSIabmRBG5FvgOeAE4BSgB3gYeAH4J3OQtxwDvAFuA7wEvAt+KyDU1DUxVn1bVu1T1VVVdUdP9w7jXW95dnsC899mEu+oDuCdMs+I9uGT9YHkC8/YrBW7ETYx5i4i0iEOMMXk3pFfiqb1a06yhPeCcCA1zszmzT/BYiqFXwcaY1Kk2iYnIdNxcXNnAH4H+qtpOVS9V1d+o6p+8hPMnVb1PVS9R1bbAAOAR3NXecyLyVQI/R3WfoSNwPHAQeC10u6p+DKwD2gInBuyXB5znrb4YZr9vgS+BPOD8uAcehqpWuS9zXn+b/DKRzg/5ft9dsNF6KRrjE7FciXUC7gC6eFdGC2M5sKouUNVfAl1w97Y61z7MOhvkLReq6r4IdWaE1AXoDTQGtkW5Ggy3X8J8s6mUb0uCx0o8y5oSE2rY0VV7KRavsgefjfGDWJJYD1X9q6oeAhCRwTXpkaeqh1T1caBHbYOMg27eclWUOqtD6gb+ezWRhdsvYUKbsk7p2drGSkywhrlVH3y2JkVj/KHaJBbmyqUY12GjRqJcASVDvrfcE6VOqbdsGof9gojITV5X/uKSkpKogVYn9H7Y+cdaU2IyVG1StAefjfGD2j4nVm3/Yu+KrXctj1+vqOpTqlqkqkUFBQW1Ps7yzbv5ZlNpxXpOlnB2P2tKTIbQB5837TrATBtL0ZiUS+TDzrcDMd0/S4Ly//mbRKlTftW1Ow77JcSk+RuD1k/q0YoWjfMS/bYG16R4Rh+bnsUYv4m1i/1qEXlFRH5eg/0iPUeWCiu9ZZcodTqF1A38d7ROKeH2S4iqDzhbU2IyXRDyQLk1KRqTerFeiXXEPej8sLd+rYisEpE3vKGXzhORwpB9euFGtvCD2d6yn4g0ilBnSEhdgCXAPqCliETqmDI0zH5x921JKUs2Vl7sZWcJ5/SzUTqS6bSQ6Vk27Tpg07MYk2KxJrFC4GLgIW+9DHcFMgr4PfBvYIOIrBGRSSLyMe55q3lxjrdWVHUNbjisPFwyDiIiw3GJeiPuua/y/Q4C73qrV4fZrztu9JKDuNE7EubdBcFNiSd2b0nLJtaUmEwNc7M545jQB583RqhtjEmGmJKYqm5R1QmqWj7qxT+BdsBI4Le40TvWAB2Ac4FTcfeTfhv3iKMQkTHeYMFjwmwuL3tQRHoG7FNIZW/LsWHGThyLG0rrbhEZGrBfPvAs7jt8QlV3xOtzhFNl2hXrlZgSodOzWJOiMalVmwGAfwiUesM1TfJeAIhIS9xIHXnAzLoMAiwigwnuyt/XWz4gIr8sL1TVEwPqtMM9oFzlf3hVHS8iT+KGmJovIlOoHAC4GfAWbiDg0P1miMg9uAGAvxCRj4AduHEcC4HpwH21/ZyxWLV1DwvX76pYzxKsKTFFTutdSKPcbPYdOgLAhp37mb1mB8d3OSrFkRmTmWqcxFR1XJRt23CjvMdDM+CEMOW9antAVb1FRD4DbsUloWzcfa9ngSfDjWDv7feQiMwDfoG7d9YQ+Bb4M/C/qnqgtjHFIrTJami3lhQ0TckMMBmvUZ5rUgwc+mvS/A2WxIxJEd9OQKWq06hh70ZVvQG4oZo6LwEv1SKe94D3arpfPFR5wNl6JabU+f3bBSWxd+dv4L7z+5Bl07MYk3S+nRTTOGu27WXe2p0V6yJwrjUlptTpxxTQMLfyV2f9zv3MWZvQW6LGmAhiGcX+ARFpXpc3EZHmIvJAXY6RqUKvwoZ0aUlhs4YpisYANM7LqdpL0R58NiYlYrkSuxs3J9h/iUiNRqIXkc4icj/u/tFdtYgv43U6qjFDu7VEvJaq82wGZ18INz2LqvVSNCbZYrkndgquA8N/Af8pItOAD3HPUy0GtqrqYW/m51a4XoQnAWfiZj3OAr7GTediaui8/u04r387Nu/az/sLN3K2NSX6whnHFNIwN4v9h1xfoHU79jFnzQ4GdbYOHsYkU7VJTFW/AoaKyFXAz4AzgNMD64jIASCwu1z5He6vgMdU9ZX4hJu5Cps15NqTuqY6DONpnJfD6b0Lgx5CnzR/gyUxY5Is5o4dqvqSqg7FDbM0Bncltg+XsBp6y73AZ8DvgMGqerIlMFNfhTYpTppvTYrGJFttnhMrxs0pBoCINAaaAztSPGeYMUl1xjGFNMjJ4sDhyibFuWt3MrBTixRHZkzmqHMXe1Xdq6obLIGZTNOkgWtSDGQzPhuTXLVOYl7Pwx+IyI9F5AIRsdkZTcY5f0Bwk+LEeRusSdGYJKpxc6KIZAGPALcQkgRFZDpuGKY34hOeMf52xjGF5OVkcTCgSXH+up0M6GhNisYkQ22uxH6Bm7U5Gzdz80RcJ4/9uOlXXhOR10QkLvOEiMhVIvKpiOwUkVIRKRaRW71kGusxuoqIxvgaFrLv/dXU3x+Pz2nSU36DHE47uiCobKI1KRqTNNVeiYlIX1VdFFD0Y+AAcIk3nmB5vca4qVl+B1yGS3KX1iU4Efkr7opvP+7ZtPJR5x8HRojI6EiD9oYoBZ6Lsr0vbmDf3cDMCHXmAnPClB+K4f1NPXbBgHZ8sGhTxfqk+Ru459xjELGxFI1JtFiaExeIyE7cA8vTgW7A+4EJDFwHD+BVEXkHNyr8lSJyZW272IvIZbgEthEYpqrLvPI2wFTgEtwV4WPVHUtVtxBlYGARKZ9O5l+quidCtbdU9f5Y4zeZY0SfNkFNimu27WPBul3071in0dqMMTGIpUluJtAYOAv4DS7xnS8iK0XkdRH5tYicIyKtAVR1P3AjsAE3d1dtlU/AeXd5AvOOvynguPfUpFkxHBHpAJzjrT5Tl2OZzJTfIIfh1qRoTEpUmwBUdQhubq+TgZ8DZbgmtM64q6E/4CbG3FSe2HD3zdYAg2sTlIh0BI4HDgKvhYnpY2Ad0BZ3H64ubsB9DwtVdXodj2UyVOiMz5PmWy9FY5Ihpt6J3qSPXwFficgvcB06fgAU4ZLN8d6/u1CZ3BRARJYBM8pfqvpZDG85yFsujPL82Qygg1f3i1g+RwQ3eMvqrsIGi8iDwFHANlzT6kRVPViH9zb1xIg+wb0UV2/by8L1uzi2gzUpGpNItZkU8wPgWqCrqk4GJpdvEJFWuIQ2BNek2B3o4b2+j7uKi+U9u3nLVVHqrA6pW2MiMhzoibvi+2c11S/0XoHWisg13pWhyWBNG+YyrFcBU1ZsIo4AAB06SURBVBZXdvCYOH+DJTFjEqw295PGAoeBd0VkZOAGVd2qqh8AD+B6BO7GXS2N8sqmxPge+d4yUicLvOMDNI3xmOH80Fu+43X+CGcF7v7cQNzwWgW4QZA/BjoCk0RkQB1iMPXEBQOCZxiwJkVjEq82YycuE5Ef47qsvy0inwJv4Jr0tuGaE+8E+gOTVXUD8I738g0RaQaM9lafjVRPVcNdoU0FporIeNzjBA/gHi+I9F43ATcBdO5coynZTBoZ0acNedlZHDzimhRXbbUmRWMSrVY9+1T1JVzT2nrcnGGP4O4RLcM9z3UxboT7X9cyrvKrrCZR6pRfre2u5Xt8H9frci3wfi2P8TtveZaI5EaqpKpPqWqRqhYVFBREqmbSXLOGuZzaq3VQmY2laExi1bp7uvecWA9cx4g3ccngCLALmACcoqqzann4ld6yS5Q6nULq1lR5U+K4GB+YDmeJt8wDWkeraDJD1elZrEnRmESqTceOCl7PvOe9VzzN9pb9RKRRhB6KQ0LqxkxE+gIn4HpQ/qN2IQJuJutypRFrmYxxZt825GYLh464xLVy614WbdhFv/bWpGhMItR5KpZEUNU1wCzcFc7lodu9XoUdcaN5fFmLt/iRt5yqqt/WNk7gCm+5VFVr26xp6pHmjXI5tVdwk7E1KRqTOL5MYp4x3vJBEelZXigihcAT3urYwKZAEblNRJaISMQrQ+/e1TXeatRnw7zpZq4SkQYh5SIi1wbE+EhMn8hkBJvx2ZjkqVNzYiKp6ngReRI3xNR8EZlC5QDAzYC3cAMBB2oN9MZdoUUyEigEduB6VUbTEngR+JuIzMJ1ZGkK9KPy+bTHVfXvsX4uU/+dFdKk+N2WPSzZuJs+7ZqlODJj6h8/X4mhqrcAV+OaFofjxjhcDtwGXKaqR2px2PIOHS954zxGswZ4GDd+ZA/c825n4b63V4ARqnp7LWIw9VjzRrl8r6f1UjQmGcSaOZKrqKhIi4uLUx2GSbDXitfwq/HzKta7FzThw58Pt+lZjKkFEZmpqkXhtvn6SsyYdHV237bkZlcmrG9L9rB0k/X9MSbeLIkZkwDNG+dySmiT4jxrUjQm3iyJGZMgob0UJ9qDz8bEnSUxYxLk7L5tyMmqbFJcUbKHbzbZM/HGxJMlMWMSpEXjPE4OaVK0GZ+NiS9LYsYk0AX9q07PYoyJH0tixiTQ2X3bkh3QpLh8cynfWC9FY+LGkpgxCXRUkzxO7tEqqGzC3PUpisaY+seSmDEJduGA9kHrb81ZZ70UjYkTS2LGJNi5/duSl1P5q7Zm2z5mrd6ewoiMqT8siRmTYM0a5nJmn8Kgsjdnr0tRNMbUL75PYt5UKJ+KyE4RKRWRYhG5VURqFLuI3C8iGuUVdTDgeMVhMtOogR2C1v89bwMHD9d2QnFjTDnfTsUCICJ/BW4B9gMfUjkVy+PACBEZHTifWIzmAnPClB9Kchwmg5zWu5AWjXPZsdf9mO3Ye4hPvinhzL5tUhyZMenNt0lMRC7DJY6NwDBVXeaVtwGmApcAtwOP1fDQb6nq/T6Iw2SQvJwsLujfjhenr64oe3POOktixtSRn5vC7vWWd5cnDgBV3YSbKBPgniQ05/klDpPmRg0KblKcsmgTu/dHbAAwxsTAl//xikhH4HjgIPBa6HZV/RhYB7QFTqzvcZj64fjOR9HxqEYV6wcOl/HegmiTkBtjquPLJAYM8pYLVXVfhDozQurGarCIPCgiT4nIWBG5RETyUhCHyTBZWVKlg8dbc6yXojF14dck1s1bropSp/zmQrcodcK5ELgL+AlwN/AGsEJEhic5DpOBRg0KfvD5ixVb2bgzasdYY0wUfk1i+d5yT5Q65XNaNI3xmCtw97cGAs2BAuAM4GOgIzBJRAYkIg4Rucnrkl9cUlISY7imPupZ2JRjOzSrWFe1YaiMqQu/JrG4U9V/qupYVZ2rqrtUdYuqTlXV04DXgcbAAwl676dUtUhViwoKChLxFiaNhDYp2oPPxtSeX5NY+dVNkyh1yq+S4jEk+O+85VkikpvCOEwGuOi49gQMbM+iDbtYsnFX6gIyJo35NYmt9JZdotTpFFK3LpZ4yzwgcBbDZMdhMkBhs4acEjJZ5isz1qQoGmPSm1+T2Gxv2U9EGkWoMySkbl0EzpUROH98suMwGeLyok5B62/OXseBw0dSFI0x6cuXSUxV1wCzcFdGl4du93oSdsSNovFlHN7yCm+5VFUrmgVTEIfJEGf3bUOLxpUt1zv2HmLyok0pjMiY9OTLJOYZ4y0fFJGe5YUiUgg84a2ODRyzUERuE5ElIvJ84IFEpLM3gG+DkHIRkWsD3uuReMRhTHUa5mZX6eBhTYrG1Jxvk5iqjgeexI2GMV9EJojIG8AyoC/wFm4A3kCtgd5A55DylsCLQImITBORl0RkAq7b/fNAI+BxVf17nOIwplpXDgluUvxs+RbWbt+bomiMSU++TWIAqnoLcDWuSW84cA6wHLgNuExVY72JsAZ4GJgJ9ABGAWfhPv8rwAhVvT0JcRhToU+7ZhzXsXnFuiq8Vrw2hREZk37EpklPrqKiIi0uLk51GMYnXpy+ivveXFCx3r55Qz69+wyyA/vgG5PhRGSmqhaF2+brKzFj6ruLjmtPo9zsivX1O/fz2fItKYzImPRiScyYFGraMJfz+7cLKntlxuoItY0xoSyJGZNioR08Ji/axNbSAymKxpj0YknMmBQb0vUoureuHNns0BG18RSNiZElMWNSTES4IuRq7JUZa7BOV8ZUz5KYMT5w6eAOQT0Sl20uZfp321IYkTHpwZKYMT5Q2LQhZ/dtE1T2/JcrUxKLMenEkpgxPnHdSV2D1t9fuIkNO/elJhhj0oQlMWN84sTuLTm6TX7F+pEy5eXp1t3emGgsiRnjEyLCtScGT1330terbYoWY6LwfRLzRp//VER2ikipiBSLyK0iEnPsIpIlIieLyB9E5AsR2S4ih0Rkk4hMEpFRUfa9X0Q0ymt/fD6pMXDJ4I7kN8ipWN9SepBJ8zekMCJj/C2n+iqpIyJ/BW4B9gMfAoeAEbhR40eIyOgYp0DpDnzu/Xsb8DWw3Ss/DzhPRMYBP9TI/ZrnAnPClB+K7dMYU738BjmMPr4j475YWVH2f598x6iBHRCx8RSNCeXbJCYil+ES2EZgmKou88rbAFOBS4DbgcdiOJwCH+FGsp8cOOq8N7HlROAG4BPgHxGO8Zaq3l+bz2JMTVx/clee+3Il5X9OLdqwiy+/3crJPVqnNC5j/MjPzYn3esu7yxMYgKpuAm72Vu+JpVlRVVeo6ghVfS902hRV/RgY661eE4e4jamTbq2bcGaf4O72T3/6XYqiMcbffJnERKQjcDxwEHgtdLuXeNbhJqo8MQ5vOdtbdozDsYyps5+c2j1o/aMlm1m+eXeKojHGv3yZxIBB3nKhqkZ6UGZGSN266OUto91BHywiD4rIUyIyVkQuEZG8OLy3MVUM6XpU0ISZYFdjxoTj1yTWzVuuilKn/AGablHqVEtEGgN3eKuvR6l6IXAX8BPgbuANYIV3T82YuBIRfhxyNfb6rLX28LMxIfyaxMqf+NwTpU6pt2xax/d6ApcIFwFPhdm+And/biDQHCgAzgA+xjU/ThKRAdHeQERu8h4NKC4pKaljuCZTnHdsWzq1bFSxfuiI8vePv01hRMb4j1+TWFKIyH8C1wM7gStUtcokTqr6T1Udq6pzVXWXqm5R1amqehruyq0x8EC091HVp1S1SFWLCgoKEvBJTH2Uk53FzcN7BpW9/PVqSnbbXGPGlPNrEiu/ymoSpU751Vqt7naLyM+B33nvdZ6qLqzFYX7nLc8SkdzaxGFMNJcd34G2zRpWrB84XMYzn9m9MWPK+TWJrfSWXaLUKZ+AaWWUOmGJyO3AH4F9wEhV/bKmx/As8ZZ5gD3EY+KuQU42Px0efG/sn1+uZNueg6kJyBif8WsSK+/y3k9EGkWoMySkbkxE5Fbgz7hRQC7yuuvXVquAf5dGrGVMHXx/SGda51d2hN1z8AhPTluewoiM8Q9fJjFVXQPMwl3hXB663esR2BE3mkfMV1Ei8h+4IasOAKNUdUodQ73CWy5VVXuIxyREo7xsfjqsR1DZc1+usp6KxuDTJOYZ4y0fFJGKu9siUojrUQgwNnDsRBG5TUSWiMjzoQcTkZ94+x0ALlHV96sLQEQ6ewMQNwgpFxG5NiDGR2rywYypqWtP6hJ0b+zg4TL+/OGyKHsYkxl8O3aiqo4XkSdxQ0zNF5EpVA4A3Ax4C3dVFag10Bt3hVZBRAYCfwcE+A64UkSuDPO2W1T1lwHrLYEXgb+JyCxgPa5Lfz8qn097XFX/XusPakwMGuZmc8eIXvz6zfkVZa8Wr+Unp3ane0F+lD2Nqd98m8QAVPUWEfkMuBUYDmTjOlM8CzwZ4wj2AC1wCQzgGO8VziogMImtwQ0aPAToCQzFXb1uBF4BnlLVj2L+QMbUweVFHXnqkxWs3LoXcJNmjnl3Cf93XVGKIzMmdSTyzCMmEYqKirS4uDjVYZg0NWHuem5/Obgv0/M/HMqwo+35Q1N/ichMVQ3715qf74kZY0Jc0L8dgzu3CCr77wkLOXQk1kYJY+oXS2LGpJGsLOH+i/oROD/mipI9PBcwiaYxmcSSmDFpZkDHFlxxfKegskcmf8Pa7XtTFJExqWNJzJg09Ktze9O0QWW/rD0Hj3DvG/Oxe9wm01gSMyYNtc5vwN3nBXey/XTZFl4rXpuiiIxJDUtixqSpq4Z25sTuLYPKfj9xkTUrmoxiScyYNJWVJTx42QAa5WZXlO3ef5jbXprNwcPWW9FkBktixqSxLq2acNe5vYPK5qzZwZh3F6coImOSy5KYMWnu+pO6MuKYwqCyf3y+kglz16coImOSx5KYMWkuK0v44xXH0aFF8KxFv3h1Ll+u2JqiqIxJDktixtQDLRrn8fhVg8jNrnwK+uCRMm56vphF63elMDJjEsv3ScybCuVTEdkpIqUiUiwit4pIrWIXkXNF5AMR2SYie0VkgYjcFzrdSpj9ThCRN0Vks4jsF5FlIvKQiDSv3SczJr4GdT6KMZcOCCrbfeAw1zwznTlrdqQoKmMSy9dJTET+ipsKpQj4FJgMHI2bgmV8TROZiNwFvAucgZt0cyJQCPwBmCYijSPs9wPgc2AU8A3wNm7Czl8Bxd4cZ8ak3OjjO1bp6LFtz0F+8NRXTF26OUVRGZM4vk1iInIZcAtu2pMBqjpSVS8BegGLgUuA22twvCJgLLAXOEVVz1TVy4HuwCfAicD/hNmvI/AMbiqXUar6PVW9EuiBm46lJ26uMmN84ebhPbjxlK5BZfsOHeFH42bw0HtLOHD4SGoCMyYBfJvEgHu95d2qWjGFrapuwk2UCXBPDa7G7sElogdVdXrA8UqBG4Ey4BYRaRGy38+ARsBzqvp2wH6HgZuAXcAoEekb8yczJoFEhN+O7Mstp/UIKi9TeGLaCi5+/HOmLtlMWZkNUWXSny/nE/OuftYAB4EWqrovTJ21QAfcVdUX1RwvD9gONAZ6quqKMHU+A04BrlbVlwLKl+Ouus5U1Q/D7PcCcDVwn6o+UN1ns/nETDI998VK7p+wkHC/5j0L8xk5oB3HdzmKjkc1JidLyA54SdVdjKmz5o1yycmu2fVTtPnE/Dqz8yBvuTBcAvPMwCWxQUDUJAb0xiWwbeESWMDxTvGO9xKAiDTDJbDy7ZH2uzogZmN84/qTu9KtdRPuGj+Pjbv2B21bvrmUR6csi7CnMYkx5efD6VmYH7fj+bU5sZu3XBWlzuqQurEcb3WUOuGO19Vb7lDVSP2UaxKHMUk37OgC3v/ZMC4d3CHVoRgTd35NYuVpek+UOqXesmkCjxeXOETkJu/RgOKSkpKogRqTCM0b5/KnKwYy4bbvcdFx7cnOssZCUz/4tTmxXlHVp4CnwN0TS3E4JoP179icP/9gEPeefwyfLtvCrFXbmbd2J3sPHuZwmXKkTCuWxiRCvP+A8msSK7+6aRKlTvlV0u4EHi/ecRjjC+2aN+KKok5cUdSp+srG+JhfmxNXessuUeqU//atjFIn9Hida3i88ntyLbxOHnWNwxhjTBz5NYnN9pb9RKRRhDpDQupGswTYB7QUkR4R6gwNPZ6q7gTKezMOqbJHhP2MMcYkhy+TmKquwQ0LlQdcHrpdRIYDHXGjeXwZw/EO4oabAtcdPvR43YGTcM+lTQzZXP6Ac7j9mgEXeqtvVheHMcaY+PJlEvOM8ZYPikjP8kJvnMInvNWxqloWsO02EVkiIs+HOd5YQIG7RWRowD75wLO47+IJVQ0dKfVR3FXc9SJyUcB+ObjhppoBb6nqolp+TmOMMbXk2ySmquOBJ4G2wHwRmSAibwDLgL7AW7iBgAO1xj3YXOXel6rOwA091Rj4whvJ/lVcc+FwYDpwX5j91gA/wiXAt0TkExH5F7Ac+L63/GndP7Exxpia8m0SA1DVW3DNeLNwieYcXNK4DbhMVWs0kqmqPgScB0zF3eO6ENgC/AYYrqp7I+z3Mm40j3eAPrjBhw8DDwNFqmrDgxtjTAr4cuzE+szGTjTGmJqJNnaiJbEkE5ESog+nFU1r3JVjukr3+CH9P0O6xw/p/xks/prroqoF4TZYEksjIlIc6a+RdJDu8UP6f4Z0jx/S/zNY/PHl63tixhhjTDSWxIwxxqQtS2Lp5alUB1BH6R4/pP9nSPf4If0/g8UfR3ZPzBhjTNqyKzFjjDFpy5KYMcaYtGVJLEVE5CoR+VREdopIqTfz860iUqtzIiLnekNpbRORvSKyQETuE5EG8Y7de7+4xC8i94uIRnntT0DsvUXkThF5wRtrs8x7r9F1PG5cz2mU94lr/CIyrppzsCTO8eeKyAgR+aP3He0SkYMisk5ExovIaXU4dsLPQSLiT/Y58N7zdhF5VUQWi8hWETkkIiUiMkVErhGRGs9eKSJZ3vdd7H3/O73z8YN4x1/Or5Ni1msi8lfgFmA/8CFwCBiBGwtyhIiMDhzYOIbj3QU8CBwBpgHbccN0/QEYKSIjIg2p5Yf4PXOBOWHKD9Ul1ghuBu6M5wET9J1EEvf4PZ/jhnULtSHO7zMcmOz9eyPwCbAHNybqZcBlIvJ7Vf1tTQ6axHOQkPg9yToHAHcDhcAC4AvcZ+gCnIH73kaLyKWxfmcikg28AVwE7AI+ABp4x3pJRE5U1fj/3KqqvZL4wv2QK+6HsldAeRtgkbftzhocrwgow/0AnhBQng987B3vER/Hf7+3z/1JPAc/Bh4CrgB64BK/AqP98J2kIP5x3v43JOn7PwMYD5waZtuVuHFJFTjdj+cgQfEn9Rx47/k9oEmY8n645KzAjTU43i+8fRYCbQLKewUc7+K4f45kfWH2qjihxd7JvC7MtuEBv4hZMR5vvLfPb8Ns6467OjsAtPBp/ElPYmFiqGsSiOt3koL4k/4faDXxPO3F80y6nIM4xO+3c/CfXjwvxVg/G9jk7TMszPbrvW1fxztWuyeWRCLSETgeN/nma6HbVfVjYB1u+pkTYzheHm5UfoAXwxzvW9ykoXnA+bUOvPL94hp/fWDfSUKUz5LeMZbKPjwHNYrfpw57ywMx1j8J1zS5VlU/CbP9NVzz7hAR6RCH+CrYPbHkGuQtF6rqvgh1ZgAdvLpfVHO83rj50bap6oooxzvFO95LNQu3injHH2iwiDwIHAVsw83vNlHdrNx+lsjvJNlOF5EBuKboTcBnwGSN3728WPXylrHeB/LbOahp/IFSfg5EpBvwH97qOzHuVn4OZoTbqKp7RWQhMNB7ratTkAEsiSVXN28ZbRT71SF1Yzne6ih1anK8WN8vXvEHutB7BVorItd4f0n7VSK/k2S7LkzZIhH5vqrOT0YAItIWuMFbfT3G3XxzDmoZf6CknwMRuRHX5JqLu3o8Gddz/QFVfTPGw8R6DgYS53NgzYnJle8t90SpU+otm6bgeKl4vxXAvbgf7uZAAe7G+ce4X6hJ3l+mfpXsc5AIc4A7cL3r8oH2wEhcj9G+wJR4NwGFIyI5wAu4n4MPVXVCjLv64hzUIX5I7Tk4BXfP6ipgmFf2n8Dva3CMlJ0DS2ImpVT1n6o6VlXnquouVd2iqlNV9TTcX7KNgQdSG2X9pqqPqupfVHWxqu5R1Q2qOhEYCnyFu9dxbxJC+RuuO/Ya4JokvF+81Tr+VJ4DVf2xqgrud60f8Ciuw9VXItI+Ee8ZT5bEkqv8L5EmUeqU/0WzOwXH89v7/c5bniUiuXE4XiIk+ztJGu9+5Bhvtc4dg6IRkceAH+G6Yo9Q1Y012D3l56CO8UeUzHOgqvtUdZGq/gqXMI/DPWMXi5SdA0tiybXSW3aJUqdTSN1Yjtc5TseL9f3iFX91ykcpyMPNJutHK71lsr6TZCs/BwlrThSRP+Ka0kpwCWBZDQ+x0lum5BzEIf7qJPwchDHOW14Y4x+QK71l0s+BJbHkKu96209EGkWoMySkbjRLgH1ASxHpEaHO0Bocrzrxjr86rQL+XRqxVmol+ztJtvJzkJDvX0QeAn4ObAXOVNVFtThMys5BnOKvTkLPQQTbcd3sc4CWMdSf5S2HhNsoIo2BY73VuJ4DS2JJpKprcCc7D7g8dLuIDMd1ZtiIe76ruuMdBN71Vq8Oc7zuuOc3DgITax145fvFNf4YXOEtl6qqL5viUvCdJFv5OQjbdbouRGQs8Cvcf5hnqeq82hwnVecgXvHHIGHnIIphuAS2A9gSQ/0vcVeiHUVkWJjtl+N6P85Q1bh1rwdsxI5kv4DRVI4e0DOgvBA3XEuV4XGA23BXXc+HOd4QKoedGhpQnk/lSA7xHHYqbvHjmkGvAhqElAtwLbDXO95PE3xOyr+niCNe4O5LLAHGxOM78Uv8uF6hI4HskPIc3DBCR7xjnxPnmP/gHXc7cHyM+/jmHMQz/lScA9yQUyOBnDDbTsH1Glbgf0O2Pe99htvC7PdLKoedKgwo7+Wdl4QMO2XPiSWZqo4XkSdxg7jOF5EpVA5U2gx4i6o3U1vjHmyucrNYVWeIyD24AYC/EJGPcH89Dcf9Ak8H7vNp/C1xI438TURmAetx3W/7UfksyeOq+vd4xQ8gIoOBJwKK+nrLB0Tkl+WFqho4ukM73GdoF3q8Wn4nfom/K/AmsM07B5txzVf9cd28y4C7VPX9OMZ/EZU/k8uB2yMMmL5EVcfG8BmSeg4SEH9XknwOgJ7AP4Ad3ntuxP3u9aDy52kirqt9oM7eZwh3j/oR3BXchcAyEfkQd/V1JtAQ+Iuqvh3Hz+DEOyvaK+a/hK7CjVi9C3cVNRO4lTBju1E5vuC0KMc7Fzey9nbcfbKFuF+0Bn6NH/eL+hAwFdcteS9uBPKVwL+AMxIU+2lePFFfIfuM88rHxeM78Uv8uD8WHsWNYrHO+/73AcuAZ4nxKqOG8d8QS/xhfl58cQ7iHX+KzkE3XO/fqbiHkPcF/O6NB0ZF2G8aUcY6xd2ius373vd45+Ez4Kp4f4byl3hvbIwxxqQd69hhjDEmbVkSM8YYk7YsiRljjElblsSMMcakLUtixhhj0pYlMWOMMWnLkpgxxpi0ZUnMGGNM2rIkZkyGE5GVIqIi0jXVsRhTU5bEjDHGpC1LYsYYY9KWJTFjMpSI3CAiSuVsvN95zYpqzYsmXdhULMZkruXAc7i5uJoArxM8e7BfZ9M2poKNYm9MhhORlbirsW6qujK10RhTM9acaIwxJm1ZEjPGGJO2LIkZY4xJW5bEjDHGpC1LYsYYY9KWJTFjzEFvaY/cmLRjScwYs85b9klpFMbUgv3lZYx5EzgNeFFEPgB2eOV3q+rWlEVlTAzsYWdjMpyIZAG/Bq4GugENvE328LPxPUtixhhj0pbdEzPGGJO2LIkZY4xJW5bEjDHGpC1LYsYYY9KWJTFjjDFpy5KYMcaYtGVJzBhjTNqyJGaMMSZtWRIzxhiTtv4/KLQQXYAV+bAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fontsize=22\n",
    "xplot = np.linspace(0,np.pi,128)\n",
    "yplot = kernel.get_values(xplot)\n",
    "plt.plot(xplot,yplot,linewidth=4)\n",
    "plt.xlabel('t',fontsize=fontsize)\n",
    "plt.ylabel('$\\phi(t)$',fontsize=fontsize)\n",
    "plt.xticks(fontsize=fontsize)\n",
    "plt.yticks(fontsize=fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Simulate a Hawkes Process\n",
    "The background intensity is $\\mu = 10$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The generated time points: \n",
      "0.011, 0.102, 0.517, 0.599, 0.678, 0.818, 0.831, 0.841, 0.883, ..., 0.011, 3.141, 3.140, 3.138 (303 points).\n"
     ]
    }
   ],
   "source": [
    "run_time = np.pi\n",
    "mu = 10\n",
    "hawkes = SimuHawkes(baseline=[mu], verbose=False, force_simulation=True, end_time=run_time, seed=1320)\n",
    "hawkes.set_kernel(0, 0, kernel)\n",
    "hawkes.simulate()\n",
    "times = np.array(hawkes.timestamps)\n",
    "print('The generated time points: ')\n",
    "for i in range(9):\n",
    "    print('{0:.3f}'.format(times[0,i]),end=', ')\n",
    "print('...',end='')\n",
    "for i in range(4):\n",
    "    print(', {0:.3f}'.format(times[0,len(times)-1-i]),end='')\n",
    "print(' ({} points).'.format(times.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Infer the Posterior Triggering Kernel by VBHP \n",
    "Search hyper-parameters over $[0.01, 10]^2$, exploit 10 inducing points and use $[0, 1.4]$ as the triggering kernel support."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/joblib/parallel.py\u001b[0m in \u001b[0;36mretrieve\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    832\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'supports_timeout'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 833\u001b[0;31m                     \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_output\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjob\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    834\u001b[0m                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/joblib/_parallel_backends.py\u001b[0m in \u001b[0;36mwrap_future_result\u001b[0;34m(future, timeout)\u001b[0m\n\u001b[1;32m    520\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 521\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mfuture\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    522\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mLokyTimeoutError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/concurrent/futures/_base.py\u001b[0m in \u001b[0;36mresult\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    426\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 427\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_condition\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    428\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/threading.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    294\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 295\u001b[0;31m                 \u001b[0mwaiter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    296\u001b[0m                 \u001b[0mgotit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: ",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-c7d132f1f4a3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;31m# Parallel\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mjoblib\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mParallel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdelayed\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 14\u001b[0;31m \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mParallel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_jobs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m7\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdelayed\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfit_vbhp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtimes\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mzmax\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrun_time\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'sin'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mga\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0malpha\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msupport\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mga\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mgammas\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0malpha\u001b[0m \u001b[0;32min\u001b[0m \u001b[0malphas\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/joblib/parallel.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, iterable)\u001b[0m\n\u001b[1;32m    928\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    929\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mretrieval_context\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 930\u001b[0;31m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mretrieve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    931\u001b[0m             \u001b[0;31m# Make sure that we get a last message telling us we are done\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    932\u001b[0m             \u001b[0melapsed_time\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_start_time\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/joblib/parallel.py\u001b[0m in \u001b[0;36mretrieve\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    853\u001b[0m                     \u001b[0;31m# scheduling.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    854\u001b[0m                     \u001b[0mensure_ready\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_managed_backend\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 855\u001b[0;31m                     \u001b[0mbackend\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabort_everything\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mensure_ready\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mensure_ready\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    856\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    857\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexception\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTransportableException\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/joblib/_parallel_backends.py\u001b[0m in \u001b[0;36mabort_everything\u001b[0;34m(self, ensure_ready)\u001b[0m\n\u001b[1;32m    536\u001b[0m         \"\"\"Shutdown the workers and restart a new one with the same parameters\n\u001b[1;32m    537\u001b[0m         \"\"\"\n\u001b[0;32m--> 538\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_workers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshutdown\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkill_workers\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    539\u001b[0m         \u001b[0mdelete_folder\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_workers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_temp_folder\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    540\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_workers\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/joblib/externals/loky/process_executor.py\u001b[0m in \u001b[0;36mshutdown\u001b[0;34m(self, wait, kill_workers)\u001b[0m\n\u001b[1;32m   1094\u001b[0m                     \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1095\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mwait\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1096\u001b[0;31m                 \u001b[0mqmt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1097\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1098\u001b[0m         \u001b[0mcq\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_call_queue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/threading.py\u001b[0m in \u001b[0;36mjoin\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m   1054\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1055\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1056\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_wait_for_tstate_lock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1057\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1058\u001b[0m             \u001b[0;31m# the behavior of a negative timeout isn't documented, but\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/threading.py\u001b[0m in \u001b[0;36m_wait_for_tstate_lock\u001b[0;34m(self, block, timeout)\u001b[0m\n\u001b[1;32m   1070\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlock\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# already determined that the C code is done\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1071\u001b[0m             \u001b[0;32massert\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_is_stopped\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1072\u001b[0;31m         \u001b[0;32melif\u001b[0m \u001b[0mlock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mblock\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1073\u001b[0m             \u001b[0mlock\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrelease\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1074\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_stop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "nparams = 8\n",
    "alphas = np.logspace(-2, 1, nparams)\n",
    "gammas = np.logspace(-2, 1, nparams)\n",
    "\n",
    "nz = 10\n",
    "zmax = 3.14\n",
    "support = 1.4\n",
    "\n",
    "# Non-parallel\n",
    "# l2s =[fit_vbhp([hp], nz, zmax, run_time, 'sin', False, ga,alpha,support) for ga in gammas for alpha in alphas]\n",
    "\n",
    "# Parallel\n",
    "from joblib import Parallel, delayed\n",
    "res = Parallel(n_jobs=7)(delayed(fit_vbhp)([times], nz, zmax, run_time, 'sin', False, ga,alpha,support) for ga in gammas for alpha in alphas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Select the posterior triggering kernel by maximizing the marginal likelihood and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mls = np.array([e[0] for e in res])\n",
    "ind = np.argwhere(mls.max()==mls)[0,0]\n",
    "paramslin = res[ind][1]\n",
    "mls = mls.reshape((nparams,nparams))\n",
    "ga_alpha = np.argwhere(mls.max()==mls)[0]\n",
    "gamma = gammas[ga_alpha[0]]\n",
    "alpha = alphas[ga_alpha[1]]\n",
    "xplot = np.linspace(0,3,128).reshape((1,-1))\n",
    "fontsize=28\n",
    "params = unlinearise_params(paramslin)\n",
    "z = np.linspace(1e-6,zmax,nz).reshape((1,-1))\n",
    "taus = [np.array([(times[0, i] - times[0, :i])[times[0, i] - times[0, :i]<=support]]) for i in range(1, times.shape[1])]\n",
    "precomp = equations.precompute(z, gamma, alpha, run_time, times[0],taus)\n",
    "mutilde, sigmatilde2, Sigmatilde = equations.predictive(params.L, params.m,\n",
    "                                   precomp,equations.precompute_predictive(xplot,z,gamma,alpha))\n",
    "ticksize = 28\n",
    "quantiles = squared_normal_quantiles(mutilde, sigmatilde2 ** 0.5, [0.10, 0.5, 0.9], double=True)\n",
    "wp.gpplot(xplot.flatten(), quantiles[:, 1], quantiles[:, 0], quantiles[:, 2],\n",
    "          fillcol='r', edgecol='r',mulabel='VBHP')\n",
    "plt.plot(xplot.flatten(), kernel.get_values(xplot.flatten()),'-.',color='g',label='Truth',linewidth=5)\n",
    "\n",
    "plt.legend(prop={'size': 26})\n",
    "plt.xticks(fontsize=ticksize)\n",
    "plt.yticks(fontsize=ticksize)\n",
    "plt.xlabel('Time Difference',fontsize=fontsize)\n",
    "plt.ylabel('Triggering Kernel Value',fontsize=fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
